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A  Unified  Formulation  of  the  Segregated  Class  of Algorithms  for  Multi  fluid  flow  at  all  speeeds 


Abstract 

In  this  paper,  the  class  of  segregated  single-fluid  all  speed  flow  algorithms  is  extended  to 
multi-fluid  flow  simulations  using  a  unified,  compact,  and  easy  to  understand  notation. 
Depending  on  the  constraint  equation  used  to  derive  the  pressure  correction  equation,  the 
extended  multi-fluid  flow  algorithms  are  shown  to  fall  under  two  categories  denoted  in  this 
work  by  the  geometric  conservation  based  algorithms  and  the  mass  conservation  based 
algorithms,  respectively.  The  differences  and  similarities  between  the  two  categories  are 
explained.  Several  techniques  developed  to  promote  and  accelerate  the  convergence  of  these 
algorithms  are  also  presented. 


Nomenclature 


coefficients  in  the  discretized  equation  for  . 

Etp  source  term  in  the  discretized  equation  for  . 
body  force  per  unit  volume  of  fluid/phase  k. 

coefficient  equals  to  1  /  . 

df  covariant  unit  vector  (i.e.  in  the  direction  of  dj. ). 

the  D  operator. 

^W[^(*)]  modified  D  operator. 

the  vector  form  of  the  D  operator, 
the  vector  form  of  the  modified  D  operator 
]  the  H  operator. 

the  HI  operator  working  on  (|)^''^  or 

HPp[(j)^''^]  the  HP  operator  working  on  or  w^^) 

//Pp[^^*^]the  modified  HP  operator  working  on  or 

HPp[u^*^  ]  the  vector  form  of  the  HP  operator. 

HPp[u^*^]the  vector  form  of  the  modified  HP  operator. 

]  the  vector  form  of  the  HI  operator, 
i  unit  vector  in  the  x-direction. 

inter-phase  momentum  transfer. 

j  unit  vector  in  the  y-direction. 

total  flux  of  across  cell  face  ‘f . 

diffusion  flux  of  across  cell  face  T . 
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convection  flux  of  across  cell  face  ‘f . 

source  of  mass  per  unit  volume. 

n^  contravariant  imit  vector  (i.e.  in  the  direction  of  ). 

P  pressure. 

heat  generated  per  unit  volume  of  fluid/phase  k. 

general  source  term  of  fluid/phase  k. 

volume  fraction  for  fluid/phase  k. 

gas  constant  for  fluid/phase  k. 

Sf  surface  vector, 

t  time. 

temperature  of  fluid/phase  k. 

interface  flux  velocity  .S y  of  fluid/phase  k. 
velocity  vector  of  fluid/phase  k. 

y(k)  velocity  components  of  fluid/phase  k. 

X,  y,z  Cartesian  coordinates. 

||a,  h||  the  maximum  of  a  and  b. 

Greek  Symbols 

density  of  fluid/phase  k. 
diffusion  eoefficient  of  fluid/phase  k. 
dissipation  term  in  energy  equation  of  fluid/phase  k. 
general  scalar  quantity  associated  with  fluid/phase  k. 


scalar  value  at  cell  face  'f . 
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Kj.  space  vector  equal  to 

A  p  the  A  operator. 

viscosity  of  fluid/phase  k. 

Q  cell  volume. 

thermal  expansion  coefficient  for  phase/fluid  k. 

5t  time  step. 

Y  scaling  factor. 

Subscripts 

e,  w, .  refers  to  the  east,  west, . . .  face  of  a  control  volume. 

E,W,..  refers  to  the  East,  West, . . .  neighbors  of  the  main  grid  point, 

f  refers  to  control  volume  face  f. 

F  refers  to  main  grid  point  F. 

P  refers  to  the  P  grid  point. 

Superscripts 

C  refers  to  convection  contribution. 

D  refers  to  diffusion  contribution. 

(k)  refers  to  fluid/phase  k. 

(k)*,(ky* refers  to  first,  second, . . .  updated  value  at  the  current  iteration. 

(k)  o  refers  to  values  of  fluid/phase  k  firom  the  previous  iteration. 

(k)'  refers  to  correction  field  of  phase/fluid  k. 

old  refers  to  values  from  the  previous  time  step. 
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SX 


x,y,z 


refers  to  SIMPLEX. 

refers  to  components  in  x,  y,  and  z  directions. 


Introduction 


Over  the  past  two  decades  important  advances  have  taken  place  in  CFD  centered  around 
increasing  numerical  accuracy  through  the  development  of  high-resolution  schemes 
[1,2,3,4,5,6,7,8,9,10,11,12,13,14],  and  improving  efficiency  through  devising  better  solution 
algorithms  [15,16,17,18,19,20,21],  better  solvers  [22,23],  and  increasing  use  of  multigrid 
techniques  [24,25,26,27].  While  high-resolution  schemes,  solvers,  and  multigrid  techniques 
can  be  applied  indiscriminately  to  the  simulation  of  single- fluid  or  multi-fluid  flows,  nearly 
all  developments  in  solution  algorithms  have  been  directed  towards  the  simulation  of 
incompressible,  compressible,  and  more  recently  all-speed  single-fluid  flows 
[25,27,28,29,30,31,32].  In  this  paper,  a  solution  algorithm  denotes  the  procedure  used  to 
solve  the  coupling  between  the  velocity,  pressure,  mass  fraction,  and  density  (for 
compressible  flow)  fields. 

For  the  solution  of  single-fluid  flow,  a  number  of  segregated  solution  algorithms  have  been 
developed  such  as  the  well-known  SIMPLE  [15,16],  the  SIMPLER  [16],  the  SIMPLEST 
[33],  the  SIMPLEC  [18],  the  SIMPLEM  [34],  the  PISO  [17],  and  the  SIMPLEX  [20] 
algorithms,  to  site  a  few.  Additionally,  several  techniques  have  been  advertised  to  improve 
the  performance,  facilitate  the  implementation,  and  extend  the  capability  of  these  algorithms. 
The  use  of  the  Momentum  Weighted  Interpolation  Method  (MWIM)  [35,36,37,38,39]  and  the 
Pressure  Weighted  Interpolation  Method  (PWIM)  [37,40,41,42,43,44]  that  have  enabled  the 
implementation  of  these  algorithms  with  a  collocated  variable  arrangement  is  an  example  of 
such  techniques.  Another  example  is  the  extension  of  the  segregated  pressure-based 
algorithms  to  simulating  compressible  and  all  speed  flows  that  span  the  entire  subsonic  to 
hypersonic  spectrum  [28,29,21,31,45,46,47,48]. 

On  the  other  hand  the  development  of  multi-fluid  solution  algorithms  has  lagged  behind  that 
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of  single-fluid  algorithms,  partly  because  of  the  much  higher  computational  cost  involved, 
and  partly  due  to  the  numerical  difficulties  that  had  to  be  first  addressed  in  the  simulation  of 
single-fluid  flow.  While  the  main  difficulty  in  simulating  single-fluid  flow  stems  from  the 
coupling  between  the  momentum  and  continuity  equations,  in  multi-fluid  flow  situations,  this 
problem  is  further  complicated  by  the  fact  that  there  are  as  many  sets  of  continuity  and 
momentum  equations  as  there  are  fluids,  that  they  are  all  coupled  together  in  various  ways 
(interchange  momentum  by  interphase  mass  and  momentum  transfer,...),  and  that  the  fluids 
share  space  (the  volume  fractions  sum  to  unity,  but  are  not  known  in  advance). 

Despite  these  complexities,  successful  segregated  pressure-based  solution  algorithms  have 
been  devised.  These  algorithms  can  be  divided  into  two  groups  based  on  the  constraint 
equation  used  in  deriving  the  pressure  or  pressure  correction  equation.  Among  these 
algorithms  is  the  IPSA  algorithm  (Inter-Phase  Slip  Algorithm)  and  its  variants  devised  by  the 
Spalding  Group  at  Imperial  College  [49,50,51,52]  and  the  Implicit  Multi-Field  [53]  based 
algorithms  (IMF)  developed  by  the  Los  Alamos  Scientific  Laboratory  (LASL)  group 
[54,55,56,57,58,59,60,61,62].  However,  in  contrast  with  the  widespread  information 
available  on  single-fluid  solution  algorithms,  much  less  information  is  available  about  multi¬ 
fluid  solution  algorithms,  a  fact  that  has  confined  their  implementation  to  a  small  community, 
slowed  their  development,  and  isolated  them  from  the  newer  developments  in  single-fluid 
flow  algorithms  (PWIM,  all  speed  flows,...). 

From  the  above,  it  is  obvious  that  the  segregated  class  of  algorithms  and  the  many  techniques 
developed  for  single-fluid  flow,  have  neither  been  fully  extended  nor  applied  to  the 
simulation  of  multi-fluid  flow.  The  main  objective  of  this  work  is  to  extend  the  applicability 
of  single-fluid  algorithms  to  multi-fluid  flow  simulations,  and  to  derive  these  algorithms 
using  a  unified,  compact,  and  easy  to  understand  notation  that  can  be  expanded 
systematically  to  yield  the  coefficients  of  the  pressure  correction  equation,  thus  facilitating 
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the  implementation  of  these  algorithms  to  a  wider  audience  within  the  CFD  community.  In 
the  process,  the  philosophies  behind  these  algorithms  in  addition  to  their  similarities  and 
differences  are  explained.  To  this  effect,  it  is  shown  that  all  multi-fluid  algorithms  can  be 
implemented  using  a  two-step  procedure,  whereby  in  the  first  step  a  multi-fluid  pressure  or 
pressure  correction  equation  is  derived  based  either  on  Volume  Conservation  (overall  volume 
fraction  equation)  [49,63,64,65]  or  Mass  Conservation  (overall  continuity  equation) 
[66,50,63,64,67,68,69,70].  Then,  in  the  second  step,  the  different  segregated  single  fluid 
algorithms  are  applied,  to  the  constraint  equation,  yielding  a  correction  equation  that  drives 
the  global  iterations  towards  convergence  in  a  manner  similar  to  the  pressure  correction 
equation  in  single-fluid  flows. 

In  what  follows,  the  governing  equations  for  multi-fluid  flows  are  presented  and  their 
discretization  outlined  so  as  to  lay  the  ground  for  the  derivation  of  the  pressure  or  pressure- 
correction  equation.  Then,  using  the  unified  notation,  the  different  multi-fluid  solution 
algorithms  are  described  and  the  framework  for  implementing  them  explained.  However,  it 
should  be  stressed  that  the  intention  of  the  paper  is  not  to  compare  the  relative  performance 
of  the  different  multi-fluid  algorithms,  this  would  be  a  work  in  progress,  rather,  the  aim  is  to 
unify  their  formulation. 

The  Governing  Equations 

In  multi-fluid  flows  the  various  eomponents  coexist  with  different  concentrations  at  different 
locations  in  the  flow  domain  and  move  with  imequal  velocities.  Thus,  the  equations 
governing  multi-fluid  flows  are  the  conservation  laws  of  mass,  momentum,  and  energy  for 
each  individual  fluid.  Moreover,  the  problem  specification  is  completed  by  the  introduction 
of  auxiliary  relations. 
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Conservation  of  mass 


The  volume  fraction  which  is  the  proportion  of  volumetric  space  occupied  by  the  k*'’  fluid 
(Q^'^VQ)  along  with  the  fluid  density,  and  velocity,  in  order  to  satisfy  the  mass- 
conservation  principle,  have  to  obey  the  differential  equation: 


dt 

Mass  sources  are  often  non-zero,  as  when  one  fluid  is  transformed  to  another  fluid. 
However,  summation  over  all  fluids  leads  to  the  “overall  mass-conservation”  equation: 


dt 


+ 


=  0 


(2) 


The  zero  on  the  right-hand  side  signifies  that  the  sum  of  mass  sources  (generation  and  loss)  is 
zero.  Thus 


(3) 

k 

The  mass  conservation  equations  can  also  be  thought  of  as  volume  fraction  equations  (i.e.  the 
equations  used  to  calculate  the  volume  fractions  occupied  by  the  various  phases  in  the  control 
volume).  In  this  case,  the  equation  can  be  rewritten  as: 


dt  ^ 

Conservation  of  momentum 


Denoting  the  velocity  of  the  k'’’  phase  by  then  the  momentum  equation  for  the  k*'’  phase 
can  be  written  as: 


v.(r'*y*>u'*>u'‘’)=  V.(r'*y*> Vu'‘>)+  r<‘'(-VP+  B'‘>)+  I'*’  (5) 
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Here  P  stands  for  the  pressure,  which  is  regarded  as  being  shared  amongst  the  fluids,  is 
the  body  force  per  unit  volume  of  phase  (k),  is  the  momentum  transfer  to  phase  (k) 
resulting  from  interaction  with  other  phases.  The  latter  term  vanishes  when  the  “overall 
momentum  equation”  is  formed.  Thus: 


k  [  dt  J  i  ^ 

The  inter-phase  momentum  transfer  accounts  for  the  interactions  among  the  fluids  and  can  be 
written  in  the  following  form 


l(*>  =  (7) 

m=  (phases  ^k) 

The  inter-phase  friction  force  causing  the  velocity  slip  is  an  internal  force  for  the  system  as  a 
whole. 


Conservation  of  energy 


Let  T^*^^  be  the  temperature  of  the  k*”  phase,  the  energy  equation  for  the  k'"  phase  is  given  by 


,th 


th  . 


3t  cp 


where 


(8) 


(9) 


frfa„<‘>V  fav'^V  fsu®  sv<‘>V' 

\A'~df*~Er)  I 

V  5z  dx  )  V  5z  dy  )  3  ^  ^ 

and  the  thermal  expansion  coefficient  of  the  k*'’  phase  which  is  equal  to  1/T^''^  for  an  ideal 


gas. 
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The  General  Multi-Fluid  Scalar  Equation 

A  review  of  the  above  differential  equations  reveals  that  they  are  similar  in  structure.  If  a 
typical  representative  variable  associated  with  phase  (k)  is  denoted  by  (|)^\  the  general 
differential  equation  may  be  written  as: 

(lo) 

where  the  expression  for  and  can  be  deduced  from  the  parent  equations.  The  four 
terms  in  the  above  equation  describe  successively  unsteadiness,  convection  (or  advection), 
diffusion,  and  generation/dissipation  effects.  In  fact  all  terms  not  explicitly  accounted  for  in 
the  first  three  terms  are  included  in  the  catch-all  source  term 

Auxiliary  Relations 

The  above  set  of  differential  equations  has  to  be  solved  in  conjunction  with  observance  of 
constraints  on  the  values  of  the  variables  represented  by  algebraic  relations.  These  auxiliary 
relations  express  physical  laws  of  various  kinds: 

a.  Geometric  Conservation  Equation 

The  equation  governing  the  volume  fractions  over  a  control  volume  can  be  written  as 
Xr^‘^  =  l  (11) 

k 

Physically,  the  geometric  conservation  equation  is  a  statement  indicating  that  the  sum  of 
volumes  occupied  by  the  different  fluids  within  a  cell  is  equal  to  the  volume  of  the  cell 
containing  the  fluids. 
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b.  Equations  of  State 

For  every  fluid  element  of  a  compressible  multi-fluid  flow,  an  auxiliary  equation  of  state 
relating  density  to  pressure  and  temperature  is  needed.  For  the  k  phase,  such  an  equation  is 
written  as: 

(12) 

In  order  to  represent  a  complete  mathematical  problem,  initial  and  boundary  conditions 
expressing  the  particular  situation  to  he  investigated  should  supplement  the  above  partial 
differential  equations  and  auxiliary  relations. 

Discretization  Procedure 


In  the  previous  sections  the  differential  equations  governing  multi-fluid  flow  phenomena 
were  presented  as  well  as  the  associated  auxiliary  relations.  The  task  now  is  to  present  the 
Finite  Volume-based  numerical  solution  algorithm  for  solving  these  equations  [16]. 


Discretization  of  the  General  Conservation  Equation 


The  general  conservation  equation  (10)  is  integrated  over  a  finite  volume  (Fig.  1)  to  yield  the 
following  expression: 


n  n 


(13) 


n  n 

where  Q  is  the  volume  of  the  control  cell.  Using  the  divergence  theorem  to  transform  the 
volume  integral  to  a  surface  integral,  one  obtains: 
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fi  s  S  n 

Replacing  the  surface  integral  over  the  control  volume  by  a  summation  of  the  flux  terms  over 
the  sides  of  the  control  volume  equation  (14)  is  transformed  to: 


^  ^  E  ■  s..] 

•  t  U  J 


nb=e,w,n,  s,Ub 


(15) 


nb=e,w,n,^t,b 

It  is  worth  noting  that  for  an  arbitrary  subdivision  of  the  volume  Q,  into  say  n  sub-volumes, 
the  conservation  equation  over  each  control  volume  can  be  integrated  and  the  global 
conservation  equation  recovered  by  adding  up  the  n  sub-volume  conservation  equations  since 
the  internal  surface  integral  appears  twice  but  with  opposite  signs  and  will  cancel. 

For  an  arbitrary  quadrilateral  (Fig.  2),  equation  (15)  can  be  rewritten  as: 


dt 


nb  =  e,w,  n,Syt,  b 


where  represents  the  total  flux  of  (j)  across  cell  face  nb  and  is  given  by 


(16) 


■  S„,  (17) 

Each  of  the  surface  fluxes  contains  a  convective  contribution,  ,  and  a  diffusive 


contribution,  hence: 


,(i)  _  ,tk)D  ,(i)C 

'*nb  —  '*nb  nb 


(18) 


The  discretization  of  the  diffusive  and  convective  fluxes  are  presented  along  an  east  face  of  a 
control  volume.  Starting  with  the  diffusive  flux,  its  value  along  the  'east'  side  is  given  by: 


^  •  S)^  =  (-r^  nSl 

=  -,<*T.<‘{(v/»),.()^)  ^  ^)<S.  -  ()«))>. 

where  y  is  a  scaling  factor  defined  as  [71]: 


(19) 
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1 


n  .d  S  .d 

e  e  e  e 


(20) 


Djand  dj.(Fig.  2)  are  the  contravariant  (surface  vector)  and  co variant  (curvilinear 
coordinate)  unit  vectors  respectively,  and  the  over  bar  denotes  a  value  obtained  by 
interpolation.  Defining  the  space  vector  as: 


/Ce  =  (21) 

the  diffusion  flux  becomes: 


de  Sg.d^  ^  ^ 


'J 


(22) 


Similar  expressions  are  obtained  for  other  faces.  The  convective  flux  for  the  'east'  side  is 
given  by: 


ji^c  ^  ^  =  c<*)^w  (23) 

where  Se  and  Ce  are  the  surface  vector  and  convection  flux  coefficient  at  cell  face  'e', 
respectively.  As  can  be  seen  from  Eq.  (23),  the  accuracy  of  the  control  volume  solution  for 
the  convective  scalar  flux  depends  on  the  proper  estimation  of  the  face  value  (t)e  as  a  function 
of  the  neighboring  node  values.  Using  some  assumed  interpolation  profile,  ([)£  can  be 
explicitly  formulated  in  terms  of  its  node  values  by  a  functional  relationship  of  the  form: 


(|)e  -  f(<t)NB,Ce)  (24) 

where  <j)NB  denotes  the  neighboring  node  values  ((|)e,  <t)w,  (!)n,  <|)s,  <1)t,  Ib,  <t>p,  <i>EE,  <t>ww,  (|)nn, 


<|>SS)  <|)bb5  etc...). 
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After  substituting  the  face  values  by  their  functional  relationship  relating  to  the  node  values 
of  (|),  Eq.  (16)  is  transformed  after  some  algebraic  manipulations  into  the  following 
discretized  equation: 

Ak)Ak)  _  V  A(k)Ak)  ,  jyik)  /-rx 

Ap  (pp  —  2^Ai^p<p^g  +  Bp  (2j) 

NB 

where  the  coefficients  Ap^  and  depend  on  the  selected  scheme  and  J^p^  is  the  source 
term  of  the  discretized  equation  .  In  compact  form,  the  above  equation  can  be  written  as 

=  Hp  [f^]=  ^ -  (26) 

Ap 

Discretization  of  the  momentum  equation 

The  discretization  procedure  of  the  momentum  equation  yields  a  discretized  equation  of  the 
form: 


NB{P)  m~{  phases  ^k) 


(  (m) 

{a/-u 


(k) 

P 


where: 


r«i 

U  =  I  V  I 

LwJ 


Tv  1 

I  2^^NB 


A^t'>  =  I  A' 


P 
{k),v 


NB 


Af'' 


^  ^NB  ^NB  ”  I  Z-f 


^NB 

(k) 

NB 


NB 


NB 


K’' 


=  1  B] 


.(k),y 


(27) 


(28) 


I  /  ,  ‘  ”  A  fS  I 

Lnb  J 

V,(P)=jv,(P)-j  j=|(VPX  I  (29) 

[v,(P).kJ 

In  the  above  equation,  the  inter-phase  term  is  written  out  explicitly  to  show  the  strong 
coupling  among  the  momentum  equations  of  the  different  fluids.  This  is  different  from  the 
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spatial  coupling  that  exists  among  the  neighboring  velocities  of  the  same  fluid.  One  way  to 
improve  the  overall  convergence  and  robustness  of  the  algorithm  is  to  re-write  the  discretized 
momentum  equations  for  the  various  phases  such  that: 


w 

Ip 


=  Ea' 


(k)  (k) 
NB^NB 


pVp 


NB 


m=  {phases  =^k) 


Where  now 


(30) 


+  24:!.*”  (31) 

m=  all  fluids 

For  later  reference,  the  value  of  before  the  addition  of  the  inter-phase  terms  will  be 

denoted  by  To  this  end,  the  discretized  form  of  the  momentum  equation  can  be 
rewritten  as: 


=  HPp  ]-  (P)  (32) 

where  the  body  force  and  inter-phase  terms  are  absorbed  in  the  source  term  within  the 

term,  or  as 


’  =  Hip  ]+  X  (33) 

m={  phases  ^k) 

where  the  body  force  and  pressure  gradient  terms  are  absorbed  in  the  B^^^  source  term  within 
the  HPp[u^*^  term.  In  equations  (32)  and  (33),  ,  and  HIp[u^*^  are  given  by. 


[Llp^^[u]  0  0  1 

=  l  0  Zy/^[v]  0  I 

L  0  0 


r  Q. 


0 


0 


0 

0 


1 


(34) 


L 


0 
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_ m={  phases  ^k) _ 


HPp[w^'‘^] 


ik)n 


NB 


X-4rv!S  +  S'*>-+Q, 

_ m  =  {phases  ^k) _ 


NB 


4ik),v 


HA 

NB 


NB  ^NB 


+<'■"+£2,  Zs' 

m~  {phases  ^k) 


{km)  (m) 

% 


(35) 


r Z4”'“S + '■  -  'i.‘’nr(v3’);  I 


2;4*>-vg+5<‘>.'-.<*>£j,(v/>y, 


Af 


“’iS  +  i*,*’’"  -  '•;*’n,(vpy. 


For  later  use,  modified  forms  of  the  HPp[u^*^  and  operators  are  defined  as; 


(36) 


HPp[u^*^]  = 


{kXw^ik) 


NB 


L  J 


= 


r  Df{u\ 
1-//P[1] 
0 

0 


\-HP[\\ 

0 


i-//p[i]J 


(37) 


Discretization  of  the  Continuity  Equation 


For  the  volume  fraction,  density,  and  velocity  fields  of  the  k*  phase  to  satisfy  the  mass- 
conservation  principle,  they  have  to  obey  the  following  differential  equation: 


(38) 
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This  equation  can  be  viewed  as  a  volume  fraction  equation  for  the  k  phase  in  which  case  it 
can  be  discretized  and  written  in  the  form 


r®  =  (39) 

or  as  a  continuity  equation  for  the  k“’  phase,  in  which  case  it  is  discretized  in  a  form  to  be 

used  for  deriving  the  pressure  correction  equation: 


V  V  )  Q_  A  = 

dt  J 

where  the  operator  represents  the  following  operation: 


(40) 


A>]=  E(<f/)  (“"t 

f^nb{F)=^e{F),  h'(F), /<  F),s(  F),t{F)MF) 

Discretization  of  the  Energy  equation 


The  discretization  of  the  energy  equation  follows  that  of  the  general  multi-fluid  scalar 
equation.  The  only  difference  is  the  one  pertaining  to  the  discretization  of  the  additional 
source  terms.  Since  a  control  volume  approach  is  followed,  the  integral  of  these  sources  over 
the  control  volume  appears  in  the  discretized  equation.  By  using  the  divergence  theorem,  the 
volume  integral  is  transformed  into  surface  integral  and  the  resultant  discretized  expressions 
evaluated  explicitly. 


Multi-Fluid  Solution  Algorithms 

Before  detailing  the  multi-fluid  segregated  solution  algorithms,  it  will  be  shown  that  these 
algorithms  may  be  grouped  under  two  categories  denoted  in  this  work  by  the  Geometric 
Conservation  Based  Algorithms  and  the  Mass  Conservation  Based  Algorithms.  To  justify  this 
classification,  attention  is  focussed  on  the  equations  and  variables  involved  in  a  multi-fluid 
flow  situation  with  n-fluids.  In  such  a  case,  there  will  be  n  momentum  equations,  n  volume 
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fraction  (or  mass  conservation)  equations,  a  geometric  conservation  equation,  and  for  the 
case  of  a  compressible  flow  an  additional  n  auxiliary  pressure-density  relations.  Moreover, 
the  variables  involved  are  the  n  veloeity  vectors,  the  n  volume  fractions,  the  pressure  field, 
and  for  a  compressible  flow  an  additional  n  unknown  density  fields.  It  is  clear  that  the  n- 
velocity  fields  are  associated  with  the  n-momentum  equations,  i.e.  the  momentum  equations 
can  be  directly  used  to  calculate  the  velocity  fields.  The  volume  fraetions  could  arguably  be 
calculated  from  the  volume  fraction  equations,  which  means  that  the  remaining  equation  i.e. 
the  geometric  conservation  equation  (the  volume  fractions  sum  to  1)  has  to  be  used  in 
deriving  the  pressure  equation,  or  equivalently  the  pressure  correction  equation.  This  results 
in  what  is  called  the  Geometric  Conservation  Based  Algorithm  (GCBA). 

Alternatively  the  equations  can  be  arranged  differently,  with  the  n  momentum  equations  used 
to  calculate  the  n  velocity  fields,  n-1  volume  fraction  (mass  conservation)  equations  used  to 
calculate  n-1  volume  fraction  fields,  and  the  last  volume  fraction  field  calculated  using  the 
geometric  conservation  equation  i.e. 

/'')  =  !-  (42) 

all  k^n 

The  remaining  volume  fraction  equation  can  be  used  to  calculate  the  pressure  field. 
However,  instead  of  using  this  last  volume  fraction  equation,  the  global  conservation 
equation  can  be  employed,  i.e.  the  sum  of  the  individual  mass  conservation  equations,  to 
derive  the  pressure  or  pressure  correction  equation.  In  this  case  the  resulting  algorithm  is 

called  the  Mass  Conservation  Based  Algorithm  (MCBA). 

With  this  classification,  attention  will  now  be  directed  towards  introducing  the  various 
available  and  possible  multi-fluid  solution  algorithms  grouped  along  the  above-mentioned 


tracks. 
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Part  I:  Mass  Conservation  Based  Algorithms 

The  sequence  of  events  in  the  Mass  Conservation  Based  Algorithm  (MCBA)  is  as  follows: 

•  Solve  the  momentum  equations  for  velocities. 

•  Solve  the  pressme  correction  equation  based  on  global  mass  conservation. 

•  Correct  velocities,  densities,  and  pressure. 

•  Solve  the  individual  mass  conservation  equations  for  volume  fractions. 

•  Solve  the  energy  equations. 

•  Return  to  the  first  step  and  repeat  until  convergence. 

The  above  steps,  along  with  some  of  the  techniques  that  were  developed  to  improve  on  them, 
are  detailed  next. 

Solving  for  Velocities 

In  this  first  step,  the  following  systems  of  momentum  equations  are  solved  to  find  based 
on  guessed  or  previously  calculated  volume  fraction  and  pressure  fields: 

=  k=\....n  (43) 

m=(^phases*k') 

From  Eq.  (43)  it  is  clear  that  the  term  couples  to  the  neighbouring  phase  (k) 

velocities  (geometric  or  spatial  coupling)  while  the  Dp^  term  couples  to 

m=(^phases^k^ 

the  velocity  of  all  other  phases  at  grid  point  P  (inter-phase  coupling).  Therefore,  the  rate  of 
convergence  of  the  iterative  solution  procedure  used  to  solve  the  above  system  will  greatly 
depend  on  its  capability  to  resolve  both  types  of  coupling.  The  spatial  coupling  represents  no 
problem  to  the  well-established  iterative  techniques  since  it  couples  velocities  of  the  same 
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phase.  The  inter-phase  coupling  is  however  problematic  since  it  relates  velocities  of  different 
phases.  An  explicit  evaluation  of  this  term  slows  the  convergence  rate  considerably  especially 
when  the  inter-fluid  momentum  transfer  terms,  represented  by  g^™"^ ,  are  large.  To  accelerate 
convergence,  the  Partial  Elimination  Algorithm  (PEA)  [68]  and  the  Simultaneous  solution  of 
Non-linearly  Coupled  Equations  (SINCE)  [67]  Technique  were  developed. 

Improvement  #1:  The  Partial  Elimination  Algorithm  and  The  Simultaneous 
solution  of  Non-linearly  Coupled  Equations  Technique 

The  central  idea  in  the  Partial  Elimination  Algorithm  (PEA),  applicable  for  two-fluid  flow,  is 
to  render  the  discretized  momentum  equations  more  implicit  by  de-coupling  the  two  sets  of 
equations.  This  is  achieved  in  a  straightforward  algebraic  manner  and  it  results  in  a 

modification  to  the  values  of  Ap  and  Bp  . 


For  the  case  of  two-fluid  flow,  the  momentum  equations  are  given  by: 

uy>=HipKn+i>‘;v''^u‘^’ 

'■  (44) 

uV^>  =  HIp[u<^>]+D<pVV;> 

It  is  clear  that  each  of  the  equations  (44)  contains  both  variables  Up^  and  u  p^  simultaneously. 
In  order  to  eliminate  the  velocity  of  one  of  the  phases  from  the  other  phase  momentum 


equation,  and  vice  versa,  these  equations  may  be  rewritten  using  the  PEA  as 


u 


(I) 

p 


u 


(2) 

P 


HIp[u<’^]+  D<;V''^HIp[u''>‘ 
HIp[u<'>]+  D<pV"^HIp[u<'> 


(45) 


In  terms  of  the  coefficients  this  is  equivalent  to 
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,(2) 


U 


(1) 


^  aX  + 


\  NB{P) 


\NB(F) 


uy>=. 


<NB(P) 


\NB(P) 


Ay>Ay)-(Q,yg<-V> 

Equation  (46)  can  be  cast  in  a  form  similar  to  equation  (25)  to  yield 


Ar 


k(2) 


.(12) 


V  (2) 


{mp)  , 


:  (1) 


a(J) 


;(2) 


=  E(ASuli)+ 


0) 


NB(P) 


^pg 


.(21) 


2]  ^  ^  ~py^  ^p^  pf* 

NB{P) 


+  (BP-rf>apVpP'^ 


+(B\?'^-r^^^QpVpP'^ 


^P  Z^y^NB^NBJ 

NBiP) 


'+ 


(2) 


(46) 


(47) 


The  use  of  PEA  renders  the  equations  more  implicit  (u^^  is  absent  from  the  equation  and 
vice  versa)  and  enhances  the  rate  of  convergence,  which  otherwise  would  have  been  slowed 
down  by  the  lagging  inter-linkage  of  the  two  equations. 

It  is  interesting  to  note  that  when  g  is  very  large  as  compared  to  the  Ap  coefficients  (i.e. 
=g»  Api)  the  and  equations  do  tend  to  a  common  value.  For  that 

purpose,  the  denominator  of  Eq.  (45)  is  written  in  terms  of  the  A^^  coefficients  and 
simplified  as  follows: 


A™A?>  -(n,g/  =  (a»> + n,g)(A<;> +n,gE  (n.g)" 

=Ai;>A®+n,g(A<;'+A?> 
Using  Eq.  (48),  the  equation  reduces  to 


(48) 


uy>: 


(2) 


+  1 


"^NBiP) 


0), 


'QpVpP 


^  /  ^NB^  NB  ^  ^ 


(2) 


<^NB{P) 


r^^^QpVpP 


P 


- 


(49) 
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In  the  limit  where  g  »  Af^,  the  equation  simplifies  to: 


r 


(1)  \NB(F) 


Up  = 


(50) 


A«  +  A^^> 

The  same  expression  appears  on  the  right-hand  side  of  the  correspondingly  reduced  equation 
for  .  The  implication  is  clear,  when  the  interface  mass  coefficients  are  high,  and 
are  everywhere  equal. 


The  use  of  the  PEA  technique  with  more  than  two  fluids  can  become  cumbersome.  The 
Simultaneous  solution  of  Non-linearly  Coupled  Equations  (SINCE)  method  [67],  is  a 
technique  similar  to  the  PEA  in  its  aim,  applicable  for  three  or  more  phases.  It  is  however 

different  than  the  PEA  in  that  the  equations’  spatial  coupling  through  the  HI p  terms  and 

inter-phase  coupling  through  the  are  accounted  for  in  two  distinct  steps. 

m=all  fluids^  k 


In  the  first  step,  the  inter-phase  coupling  is  resolved  by  solving  the  momentum  equations  on 

u 

terms  to  the  right  hand  side  and  treating  them  as  source  terms.  The  equation  for  a  k-fluid  flow 
can  be  re-arranged  and  written  in  the  following  form: 


0(1)  ^  DO) 

Up  i^p 

^gO».)0(«) 

+  HI^[u<'>] 

m= 

:  all  fluids 

uf  =  D<,2> 

+  HIp[u<2> 

m 

=all  fluids^! 

u(;>  = 

m 

=all  fluids^  k 

(51) 


These  equations  are  solved  using  any  desired  technique  (implicit  or  explicit)  for  the 

U  p  which  are  then 


used  to  calculate  the  inter-phase  friction  terms 
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jjW  for  all  k  (52) 

m=all  fluids 

In  the  second  step,  the  inter-phase  friction  terms  are  absorbed  into  the  standard  source  terms 
of  equation  (25)  and  the  full-field  solution  of  the  momentum  equations  is  accomplished  in  a 
sequential  manner  using  the  normal  calculation  method  via  the  following  system  of 
equations: 

ur  =  HIp[u^*'‘]+4*'  for  all  k  (53) 

m  =all  fluids 

Solving  for  Pressure  Correction 

Before  convergence,  the  velocities  calculated  from  the  momentum  equations  do  not 
necessarily  satisfy  the  mass  conservation  equations.  In  the  segregated  approach,  the  burden  of 
restoring  balance  rests  on  the  pressure  correction  equation,  which  in  this  case  is  derived  from 
the  overall  mass  conservation  equation.  Therefore,  the  first  step  in  developing  a  segregated 
solution  algorithm  is  to  derive  such  an  equation.  The  segregated  MCBA  can  be  viewed  as  an 
extension  to  the  SIMPLE  algorithm  and  its  variants,  in  which,  the  pressure  equation 
derivation  follow  the  same  pattern  as  in  SIMPLE  (or  any  of  its  variants),  and  the  corrections 
are  applied  only  to  the  velocity,  pressure,  and  density  fields.  No  correction  is  applied  to  the 
volume  fractions,  rather,  they  are  obtained  by  solving  the  individual  continuity  and  geometric 
constraint  equations. 

To  derive  a  pressure  or  pressure-correction  equation,  the  various  continuity  equations  are  first 
added  to  yield  the  total  mass  conservation  equation  given  by: 


The  derivation  starts  by  noticing  that  in  the  predictor  stage  (the  previous  step)  a  guessed  or 
estimated  pressure  field  from  the  previous  iteration  denoted  by  P° ,  is  substituted  into  the 
momentum  equations.  The  resulting  velocity  field  denoted  by  ,  which  now  satisfies  the 
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momentum  equations,  in  general,  will  not  satisfy  the  continuity  equation.  Thus,  a  correction 
is  needed  in  order  to  yield  velocity  and  pressure  fields  that  would  satisfy  both  equations. 
Denoting  the  corrections  for  pressure,  velocity,  and  density  by  F' ,  ,  and 

respectively,  the  corrected  fields  are  written  as: 

'P  =  F  +  F 

\  (55) 

pW  ^  pfkY  ^  pfkY 

Thus,  before  the  pressure  field  is  known,  the  velocity  obtained  from  the  solution  of  the 
momentum  equations  is  actually  rather  than  Hence  the  equations  solved  in  the 

predictor  stage  are: 

nf*  =  (56) 

While  the  final  solution  satisfies 


(57) 

Subtracting  the  two  sets  of  equation  (57)  and  (56)  from  each  other  yields  the  following 
equation  involving  the  correction  terms: 

uf  =  HPp[u‘'‘’'  ]  -  DfVpF  (58) 

The  velocity  and  density  fields  are  corrected  to  satisfy  mass  conservation.  Therefore,  the 

new  density  and  velocity  fields,  and  will  satisfy  the  overall  mass  conservation 

equation  if: 


X r'  P' )  )  n+A,p‘V‘v‘’.s]i- = 0 

*  i  ■*  j 

Linearizing  the  (p®u^^)  term,  one  gets 


Old 


1 


(59) 


Substitution  of  equation  (60)  into  equation  (59)  gives 


(60) 
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{k)o  (k)\  {ik)  iky 


Mrp  Pp  )-{r,  Pp  )  | 

Rearranging,  the  following  equation  is  obtained: 


=  -  A,[r<*>-(p“’V*>‘  +p<‘>u<*'  )S  j 


Using  equation  (58),  the  above  equation  becomes: 


^  {§  ^  A^  ]-  D'''Vp')s  j- 

f  (it)*!  (^)*  /  (k)  (k)^^^  1 

*  I  J 

Finally,  substituting  density  correction  by  pressure  correction,  as  obtained  from  the  equation 
of  state,  the  final  form  of  the  pressure-correction  equation  is  written  as: 


2{| 


=  -l< 


''  V'  f  n + A, [/«>“>•{/<*'•  ] 


The  second  order  correction  term  d'"'’  is  usually  neglected.  This  does  not  affect  neither 
the  convergence  rate  (i.e.  it  is  considerably  smaller  than  other  terms)  nor  the  final  solution. 


since  at  the  state  of  convergence  the  correction  fields  vanish.  For  this  reason,  it  is  neglected 


in  all  subsequent  derivations.  Furthermore,  if  the  HP|u^'‘^  term  in  the  above  equation  is 

retained,  there  will  result  a  pressure  correction  equation  relating  the  pressure  correction  value 
at  a  point  to  all  values  in  the  domain.  Eventhough  such  an  equation  ensures  that  the  corrected 
fields  will  satisfy  the  continuity  and  momentum  equations,  is  undesirable  because  it  becomes 
intractable.  To  facilitate  implementation  and  reduce  cost,  simplifying  assumptions  related  to 
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this  term  have  been  introduced.  Depending  on  these  assumptions,  different  algorithms  are 
obtained.  These  algorithms  were  originally  developed  for  single-fluid  flow  and  most  of  them 
have  not  yet  been  extended  to  multi-fluid  flow.  It  is  an  objective  of  this  work  to  perform  this 
extension. 

Extending  the  Segregated  class  of  single-fluid  flow  algorithms  to 
multi-fluid  flow  situations 

Moukalled  and  Darwish  [21]  have  recently  imified  the  formulation  of  the  segregated  class  of 
algorithms  for  single-fluid  flow.  In  their  work,  the  SIMPLE  [15,16],  SIMPLER  [16], 
SIMPLEC  [18],  SIMPLEST  [33],  SIMPLEX  [20],  SIMPLEM  [34],  PISO  [17],  and  PRIME 
[72]  algorithms  were  presented  in  details  and  the  philosophies  behind  them  explained. 
Therefore,  it  is  sufficient  here  to  present  the  symbolic  form  of  the  multi-fluid  version  of  these 
algorithms  followed  by  a  brief  discussion  of  their  merits  and  the  philosophies  behind  their 
development. 

In  the  derivations  to  follow,  the  superscripts  “old”  and  “o”  denote  values  from  the  previous 
time  step  and  values  from  the  previous  iteration,  respectively.  Moreover,  the  superscripts  *, 
***^  ****  represent  the  first,  second,  third,  and  fourth  updated  values  at  the  current 

iteration,  respectively. 


The  MCBA  following  SIMPLE  (MCBA-SIMPLE):  Symbolic  Form 

Predictor: 

=  HPp[u^''^*  ]  -  (65) 

Corrector: 


(66) 
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\uf'  =  HP^[uW']- 


Condition: 


^  ^  Kp  Pp  )  Pf..±_  Q  +  .s]|  =  0 


r  (A)»  ,  Aky 
[fp  +Pp 

c 

~!<e 

si^ 
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]  ■* 
+  A,p"( 

.afi. 

ho 

2 1-  r-‘^'d‘'Pr  +  A,p‘V“’’(''“’'D®VP')sl- 

k  ^St 


Sk)o  iky 


sjDld 


.  ,  .  .  /  (^)  (*)v 

=  -YJx  ^  ^  J 

*  [+ App*^>^*^°(HP[u^^'])s]+ 


Approximation: 


Neglect:  HPp*^’  , 

^  uf  = 

Approximate  Equation: 


+  A^  ]-  A^  p*)” 


::::> 


.  1  <* 


A  Global  MCBA-SIMPLE  Iteration 

•  Solve  implicitly  for  u'‘‘^  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D*‘‘*  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u^‘‘\  P  and 

•  Solve  implicitly  for  F’. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 
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The  MCB A  following  SIMPLEC  (MCBA-SIMPLEC):  Symbolic  Form 


Predictor: 

nf*  =  HPp[u^^*  ]  -  (74) 

Corrector: 

(u^*>**  =  =  F  +P\p^'‘^  =  (75) 

=  HPp[u‘^^“  ]  -  =  HPp[u^'‘^’  +  +  P')  (76) 

=  HP;,[u^'^^'  ]  -  (77) 

Subtracting  HPp  [l]up  ^  from  both  sides,  one  gets 

nf  -  HPp[l]u^*^  =  HP/,[u^*^']  -HP^[l]up  (78) 

/.(l-HPp[l])u^;>'  (79) 

^  HPp[u<->;-u^f]  _ 

;.<{  l-HPp[l]  l-HPp[l]  (80) 

=cl*y 

Condition: 
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Approximation: 


(83) 


Neglect:  HPju'*^ 


Approximate  Equation: 


(84) 


X  {|  rf^°(^PP'p  +  p*>“C5/)t/<*>*P'  ]-  A^  VP' 


*  I 
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•  Q+A^p)y*>“t/w*]| 


A  Global  MCBA-SIMPLEC  Iteration 


•  Solve  implicitly  for  u(k),  using  the  old  pressure  and  desity  fields. 

•  Calcuate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u®  P  and  p®. 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  MCBA  following  PRIME  (MCBA-PRIME):  Symbolic  Form 

Predictor: 


=  HPp[u®°]  -r‘*^°DV*VpP° 

Corrector: 


=HP,[u'^>**]-r<*^’D^*V,P*  =  HP,[u^‘'^’  +P') 


(86) 


(87) 


(88) 


F 
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Condition: 


*  I 
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+Ap  P*^>^*^°(hp[u^*^*  -  u<*^]+  HP[u^*^'])s]f  App*V^*^ 


Approximation: 


Neglect:  HP[u^*^*  ,  HP[u^'‘^'  , 

Approximate  Equation: 


^  M)o(-\k)  O'  .  A 
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(93) 


A  Global  MCBA-PRIME  Iteration 

•  Solve  explicitly  for  u^'‘^  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D^'‘^  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u®,  P  and  p®. 

•  Solve  implicitly  for  /'‘I 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 
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The  MCBA  following  SIMPLEST  (MCBA-SIMPLEST):  Symbolic  Form 


Predictor: 


Corrector: 


nf'*  =  HP/  [u^'^’*]+  HP/[u<*^**  ]-  r(*>'D/VpP* 

=  HP/[u<*)*  +  u<*>']+HP/[u<*'‘  +  P' 

=  HP/ [u<'^^*]+  HP/[u<*>^  ]+  HP/[u^'^* ]+  HP/[u<*^  ] 


nf  =  HPp[u<*n+  HP/[u<*>* 


pW  ^  ^(k)p, 


Condition: 
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(100) 


+ A^ P*^><^^°(hP^[u<^^*  - HP[u<*^'])s]+  A^p*>y*^'u'*^'.s| 


ADProximation: 


Neglect:  HP^[u<'^‘ -  ,  HP[u<'>'  , 


nf  =  -r<'>°D/V^F 


(101) 
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Approximate  Equation: 


“I 


A  Global  MCBA-SIMPLEST  Iteration 


J 


•  Solve  for  u®,  treating  implicitly  and  HP^  explicitly. 

•  Calculate  the  D^'‘^  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P  and 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  MCBA  following  SIMPLER  (MCBA-SIMPLER):  Symbolic  Form 
First  Predictor: 

No  predictor  stage.  Only  coefficients  of  the  momentum  equations  are  calculated. 
First  Corrector: 


Condition: 


*  I 


{(FkyAk)*\_(j.(k)  (k)>P'‘‘  1 
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Approximation: 


Neglect:  HP[u'*^'  , 

=>  nf  = 

Approximate  Equation: 
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a’ 


.(k) 
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*  I 

Apply  correction  to  pressure  and  density  fields  only. 

Second  Predictor: 

=  HPp[u^''^*  ]  - 

Second  Corrector: 
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(111) 


(112) 

(113) 

(114) 
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(117) 
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Approximation: 

Neglect: 


uW-  =  -r(*>D<*V^P" 


(119) 


Approximate  Equation: 
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(120) 
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Do  not  correct  pressure. 

A  Global  MCBA-SIMPLER  Iteration _ 

•  Calculate  the  D^‘‘^  fields. 

•  Solve  the  pressure  correction  equation  and  update  the  pressure  and  density  fieids. 

•  Solve  implicitly  for  u^‘‘^  using  the  new  pressure  and  density  fields 

•  Solve  the  pressure  correction  equation  using  the  new  velocity  fields  to  obtain  a  new 
pressure  correction  field. 

•  Correct  u®. 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 
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The  MCBA  following  PISO  (MCBA-PISO):  Symbolic  Form 

First  Predictor: 

uf*  =  HPp[u®‘]  - 
First  Corrector: 

(u'*’"  =  u<‘>'  +  u“>>'  =  F  +P',p<*>‘  =  p<‘>-  +/>' 

[p“>'=CfP 
Condition: 
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Approximation: 


Neglect:  HP[u^'^^  , 
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Approximate  Equation: 
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Approximate  Equation: 
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A  Global  MCBA-PISO  Iteration 


(139) 


Solve  implicitly  for  u'  using  the  old  pressure  and  density  fields. 
Calculate  the  D^'‘^  fields. 

Solve  the  pressure  correction  equation. 

Correct  P,  and  . 

Solve  Implicitly  for 

Solve  implicitly  the  energy  equation  and  update  the  density  fields. 
Solve  the  momentum  equations  explicitly  and  calculate  the  D®  fields. 
Solve  the  pressure  correction  equation. 

Correct  P ,  and 

Return  to  step  one  and  iterate  until  convergence _ 


The  MCBA  following  SIMPLEX  (MCBA-SIMPLEX):  Symbolic  Form 


Predictor: 


=  HPp[u^''^‘]  -  pP° 

Corrector: 


(140) 


(u<*^**  =  =  F  +P',yO<*^*  = 

jnf  =  HP^[u(*>']- 


(141) 


(142) 


Condition: 
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Approximation: 
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Neglect:  and  let 

=  HPp[u^*>'  ]-  (147) 

(148) 

Assume  that  the  pressure  difference  local  to  the  velocity  is  representative  of  all  pressure 
differences i.e.  HP^[-/'^°D^*^-^V^P']  =-(VpF)HPp[/*^°D^*^^-^],  thus: 


Approximate  Equation: 
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A  Global  MCBA-SIMPLEX  Iteration 

•  Solve  implicitly  for  using  the  old  pressure  and  density  fields. 

•  Calculate  the  D®  fields. 

•  Solve  implicitly  for  the  fields. 

•  Solve  the  pressure  correction  equation  using  these  D  fields. 

•  Correct  u^'‘\  P  and  p^'‘( 

•  Solve  implicitly  for  F\ 

•  Solve  Implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


A  Unified  Formulation  of  the  Segregated  Class  of Algorithms  for  Multi-fluid  Fluid  Flow  at  All  Speeds 


41 


The  MCBA  following  SIMPLEM  (MCBA-SIMPLEM):  Symbolic  Form 

First  Predictor: 


No  predictor  stage.  Only  coefficients  of  the  momentum  equations  are  calculated. 
First  Corrector: 


=  HPp[u^*^*  ]  -  r^*^°DlfVpP’ 
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Approximation: 


Neglect:  HP[u<*^'  , 
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Approximate  Equation: 


I 


J 


Second  Predictor: 


u. 


(160) 


Second  Corrector: 


No  corrector  stage. 

A  Global  MCBA-SIMPLEM  Iteration 


•  Calculate  the  0'“^  fields  based  on  values  from  the  previous  iteration. 

•  Soive  the  pressure  correction  equation. 

•  Correct  u®,  P  and  p®. 

•  Caiculate  new  and  fieids. 

•  Soive  impiicitiy  for  u  using  the  new  fieids. 

•  Solve  implicitly  for  /‘'I 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ _ 

Discussion  of  the  MCBA  following  the  SIMPLE-like  segregated  Approach 


From  the  above  derivations  it  is  clear  that  the  main  difference  among  the  various  MCBA  is  in 
the  approximation  made  to  the  HP[u^^^  ]  term.  Since  at  the  state  of  convergence  the  pressure 
and  velocity  correction  fields  are  zero,  the  approximation  introduced  does  not  have  any  effect 
on  the  final  solution.  Rather,  it  only  affects  the  convergence  behavior. 

As  mentioned  before,  the  majority  of  SIMPLE-like  algorithms  have  not  been  extended  and 
tested  in  multi-fluid  flow  situations.  As  such,  the  only  way  to  discuss  their  merit  is  by 
assuming  that  they  behave  the  same  way  as  in  single-fluid  flow  situations. 

In  the  MCBA-SIMPLE  algorithm,  the  HP[u^*^  ]  terms  are  dropped.  Because  of  this,  the 
predicted  pressure  correction  field  is  overestimated  and  the  corrected  velocity  field  does  not 
satisfy  the  momentum  equations.  Therefore,  in  order  to  avoid  divergence,  the  pressure  field  is 
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under-relaxed.  This  under-relaxation  procedure  has  resulted  in  the  rate  of  convergence  to  be 
greatly  dependent  on  the  proper  choice  of  the  under-relaxation  factors  for  the  velocity 
components  and  the  pressure. 

The  SIMPLEC  algorithm  of  Van  Doormaal  and  Raithby  [18]  was  developed  with  the 
intention  of  alleviating  the  aforementioned  problem  through  a  better  approximation  to  the 
HP[u^*^  ]  (i.e.  neglecting  HP[u^*^  ]  rather  than  HP[u^*^  ]),  so  as  to  eliminate  the  need 

for  under-relaxing  the  pressure  field.  Consequently,  a  higher  rate  of  convergence  is  expected 
intheMCBA-SIMPLEC. 

In  the  MCBA-PRIME  algorithm  [72],  the  momentum  equations  are  solved  explicitly.  This 
explicit  treatment  is  justified  by  the  small  contribution  to  the  convergence  of  the  entire  flow 
field  by  the  iterative  sweeps  within  each  momentum  equation.  On  the  other  hand,  finding  the 
correct  solution  for  the  pressure  field  represents  the  most  important  factor  in  the  overall 
convergence.  Moreover,  the  terms  neglected  in  MCBA-PRIME 

^Ppju^'‘^  J+HP^[u®* can  become  smaller  than  the  term  neglected  in  MCBA- 

SIMPLE  ^Pp[u^’‘’  ]  if  HPp[u®  and  HPp[u^'‘^*  are  of  opposite  signs.  It  is  worth 

noting  that  HPpjii^''^  =  HPp[u^'‘^** -u®*  is  a  correction  to  satisfy  continuity,  while 

HPp[u^''^*  -  is  a  correction  to  satisfy  momentum.  Usually  the  correction  added  to  satisfy 
momentum  is  opposite  to  that  added  to  satisfy  continuity  and  hence,  the  neglected  term 

(HPp [u®  ]-!-  HPp  [u^^*  -  u®°  ]  is  smaller. 

In  MCBA-SIMPLEST,  the  coefficients  of  the  momentum  equations  are  separated  into  their 
diffusion  and  convection  parts  as: 

HPp[u<^^]=  HPf  [u<*^]+HP^[u^*>  (161) 

The  momentum  equations  are  then  solved  treating  the  diffusion  terms  implicitly  and  the 

convection  terms  explicitly.  This  procedure  is  justified  by  the  fact  that  disturbances  in  a  pure 
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diffusion  situation  propagate  instantaneously  in  all  directions,  but  their  amplitude  decays 
rapidly.  This  is  equivalent  to  propagation  of  errors  throughout  the  entire  solution  domain,  in  a 
single  iteration  by  implicit  solution  methods.  On  the  other  hand,  disturbances  in  a  pure 
convection  situation  are  propagated,  in  the  flow  direction,  at  a  finite  speed  without  any 
change  in  their  magnitude.  This  is  similar  to  propagation  of  error,  from  a  particular  point  to 
the  neighboring  grid  points,  in  a  single  iteration  of  explicit  iterative  methods.  This  justifies 
the  implicit  and  explicit  treatments  of  the  diffusion  and  convection  parts,  respectively. 
Moreover,  similar  to  MCBA-PRIME,  the  terms  neglected  in  MCBA-SIMPLEST  can  become 
smaller  than  the  term  neglected  in  MCBA-SIMPLE.  Therefore,  the  MCBA-SIMPLEST  can 
be  seen  to  be  a  compromise  between  MCBA-SIMPLE  and  MCBA-PRIME. 

As  presented,  the  MCBA-SIMPLER  can  be  viewed  as  a  combination  of  an  approximate  form 
of  MCBA-PRIME  to  compute  the  pressure  and  MCBA-SIMPLE  to  compute  the  velocities, 
thereby  combining  the  best  features  of  both  methods. 

The  MCBA-PISO  consists  of  a  predictor  step  and  one  or  more  corrector  steps.  The  use  of  the 
additional  corrector  steps  bring  the  velocity  and  pressure  fields  closer  to  satisfying  both  the 
momentum  and  continuity  equations.  Moreover,  by  following  the  sequence  of  events,  it  can 
be  easily  seen  that  MCBA-PISO  may  be  considered  to  be  a  combination  of  one  MCBA- 
SIMPLE  step  and  one  MCBA-PRIME  step,  hence  combining  the  implicitness  of  the  SIMPLE 
algorithm  with  the  stability  of  the  PRIME  algorithm. 

In  all  MCBA-SIMPLE-based  methods,  no  care  is  taken  to  ensure  that  the  rate  of  convergence 
will  not  degrade  with  grid  refinement.  This  concern  is  addressed  in  MCBA-SIMPLEX  by 
using  extrapolation  to  express  all  pressure  differences  in  the  domain  in  terms  of  the  pressure 
difference  local  to  the  velocity.  The  idea  is  based  on  the  fact  that  the  spatial  distribution  of 
pressure  difference  influence  changes  little  with  grid  refinement.  Therefore,  if  the  pressure 
difference  influence  were  restricted  to  a  control  volume,  it  would  be  appropriate  to  assume 
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that,  by  extrapolation,  the  pressure  difference  at  the  main  grid  point  adequately  represents  the 
pressure  differences  at  the  control  volume  faces.  Moreover,  the  pressure  correction  equation 
in  MCBA-SIMPLEX  involves  the  fields,  rather  than  the  fields,  which  have  to  be 
computed  by  solving  an  additional  system  of  equations  (Eq.  (150)). 

The  pressure  field  in  MCBA-SIMPLEM  is  computed  using  the  old  velocity  fields.  This  is 
nearly  equivalent  to  MCBA-PRIME,  which  does  a  good  job  in  correcting  the  pressure  field. 
By  so  doing  however,  the  velocity  corrections  will  be  at  a  disadvantage.  Therefore,  in 
MCBA-SIMPLEM  the  disadvantages  and  advantages  of  MCBA-SIMPLE  are  interchanged. 

The  Expanded  Form  of  the  Pressure-Correction  Equation 


It  is  obvious  by  now  that  the  various  simplified  pressure-correction  equations  are  similar  and 
may  be  written  as: 


^  {1  A/,  D<*>VP')S 


\Old 


1 


^  Pp  Pp  )  Q + 


(162) 


Where,  depending  on  the  algorithm  used,  and  represent  values  from  the  previous 
iteration  or  fi'om  a  previous  corrector  step.  When  discretizing  this  equation,  careful  attention 
should  be  paid  to  the  second  term  on  the  left  hand  side  that  is  similar  to  a  convection  term 
and  for  which  any  convective  scheme  may  be  used.  Adopting  the  UPWIND  scheme  [16]  the 

discretized  form  of  the  convection-like  term  A^  is: 
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+  )  +  (rflC'y‘>P  )  +  (r"’ 

+ (r'‘>'c«>)_|]^!‘',o||F;  -||-c^*',0||p;]+  (r'‘>'ci*')  o||f;  -||-c//’,o||/?] 

+  - |-i7;*’,o|p;]+  (r<‘>v;’)||c4‘’,0||/?  -||-C/'‘>,0||/s] 

Rearranging,  one  obtains: 

r  (,<•>•<?;>  )Jc^*),  oil + (r«  V/>)J|c/«>,oI  1 

l_H-(r<«-C<«)||C/;*>,0|+(r«>-C<*')||yJ*'.0|J 

(164) 

-(r<‘>-cJ*>)j|-£<‘',o||F;-(.<«-d;>;|||-c/i‘'.o||F; 


The  term  A,p’>">'(-fl?’VF)s  is  discretized  following  the  same  procedure  that  was 
used  in  discretizing  the  diffusion  flux.  Its  final  form  is  given  by: 

Wp')s)i  +  D<*Vf)s)  (165) 


Or  in  expanded  form  as, 
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(166) 


where  the  underlined  terms  account  for  the  non-orthogonal  factors.  They  are  usually 


neglected  since  their  contribution  is  small  in  comparison  with  other  terms  and  vanish  when 
the  grid  is  orthogonal.  However,  they  could  be  accounted  for  by  moving  them  to  the  right 
hand  side,  adding  them  to  the  source  term,  and  modifying  the  solver  to  explicitly  update  their 
values  after  a  solver  (not  global)  iteration.  Neglecting  these  terms,  Eq.  (166)  becomes: 
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Moreover,  the  discretized  form  ofAp[r^*^  is  given  by: 


(167) 


+  yrp\^^uy+ry°pf^uy 


(168) 


Substituting  the  various  terms  in  Eq.  (162)  by  their  equivalent  expressions  as  derived  above, 
the  final  form  of  the  pressure-correction  equation  is  written  as: 


<;>  =  4>;  +  A^F,  +  <>;  +  A[p; + a^p;  +  +  s; 


(169) 


where 
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The  corrections  are  then  applied  to  the  velocity,  pressure,  and  density  fields  using  the 
following  equations: 


P  =P° +P’  (171) 

Improvement  #2:  Weighted  Pressure  Correction 

Numerical  experiments  [73]  using  the  above  approach  to  simulate  air-water  flows  have 
shown  poor  conservation  of  the  lighter  fluid.  To  understand  this  behavior,  the  residual  error 
of  the  k*  phase  continuity  equation,  after  any  global  iteration,  that  arises  because  the 
velocity,  density,  and  volume  fi'action  fields  do  not  yet  satisfy  the  continuity  equation  is 
denoted  by  RESC^\  The  pressure  correction  equation  being  derived  from  the  global 
conservation  equation,  it  is  intended  to  correct  the  velocity  fields  so  as  to  drive  the  global 
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residual  error,  which  is  equal  to  the  sum  of  the  local  residuals,  to  zero  i.e.  RESC^’^+  RESC*^^^ 
+  ..  +  RESC^"^^0. 

In  the  presence  of  a  relatively  very  high  density  fluid  such  that  p(n)»p(k)  for  kn,  the 
residual  error  of  the  n*  fluid  will  be  of  a  magnitude  commensurate  with  the  respective  phase 
density,  i.e.  RESC^"^  is  expected  to  be  much  larger  than  RESC^*'^  for  k  n .  In  this  case,  only 
the  residual  of  the  high  density  fluid  will  be  significant  while  that  of  the  low  density  fluid  will 
be  relatively  negligible,  and  hence  the  pressure  correction  will  tend  to  drive  the  high  density 
fluid  to  conservation. 

This  problem  can  be  considerably  alleviated  by  normalizing  the  individual  continuity 
equations,  and  hence  the  global  mass  conservation  equation,  by  means  of  a  weighting  factor 
such  as  a  reference  density  (which  is  fluid  dependent)  to  give  a  volumetric  conservation 
equation  of  the  form: 


+  ^  |}-  =  0 


(172) 


In  this  case  the  pressure  correction  equation  becomes 


Solving  for  volume  fractions 

The  volume  fractions  can  be  obtained  by  solving  the  fluid  continuity  equations  (Eq.  (4)). 
However,  if  all  the  volume  fractions  are  obtained  by  solving  the  continuity  equations,  the 
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geometric  constraint  will  not  be  enforced  unless  the  appropriate  velocity  field  is  available, 
which  is  not  the  case  until  convergence.  One  remedy  is  to  solve  for  n-1  volume  fractions 
using  the  fluid  mass  conservation  equations,  and  use  the  geometric  conservation  equation  to 
find  the  last  volume  fraction  field.  Thus,  for  fluids  k=l  to  n-1,  solve 


=  (174) 

and  for  k=n  use 

=  (175) 

fluids^  n 

A  drawback  of  this  procedure  is  that  the  volume  fraction  field  of  any  phase  is  not  influenced 
by  the  volume  fraction  fields  of  other  phases  except  when  calculating  the  n*  phase.  This  can 
affect  the  convergence  rate  negatively.  Better  solutions  are  presented  below. 

Improvement  #3:  Mutual  influence  of  Volume  Fractions 

An  improvement  on  the  above  procedure,  is  to  solve  all  continuity  equations  to  obtain  all 
volume  fractions  and  then  enforce  the  geometric  conservation  constraint  on  the  resolved 
volume  fractions  using  the  following  equation: 


m-\.M  m  =  \..n 


(176) 


where  is  the  volume  fraction  resulting  from  the  solution  of  the  volume  fraction  equation 
for  fluid  k.  The  summation  is  carried  out  for  all  phases,  and  r^''^  is  the  value  of  the  volume 
fraction  for  fluid  (k)  which  is  carried  into  subsequent  calculations.  These  r^'‘^  of  course,  do 
sum  to  unity,  and  the  solution  of  each  of  the  volume  fraction  equations  affect  that  of  the 
remaining  phases  through  the  enforcement  of  the  geometric  constraint  over  all  phase  volume 
fractions  rather  than  over  the  n*’’  phase  only. 
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Improvement  #4:  Implicit  Volume  Fraction  Equations 

The  solution  of  the  volume  fraction  equations  can  be  improved  by  implicitly  accounting  for 
the  influence  of  the  volume  fractions  of  the  different  phases  on  each  other.  The  details  of  the 
procedure  will  be  presented  for  a  two-fluid  flow  and  then  generalized  to  n-fluid  flow 
situations.  For  that  purpose,  the  following  simplified  form  of  the  volume  fraction  equation  is 
considered: 

rr =»[/*'■]= 

For  the  case  of  a  two-fluid  flow,  the  volume  fraction  equations  can  be  written  as 
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The  sum  to  1  rule 


=  1  can  be  written  in  the  following  form 
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Based  on  this  equation,  one  can  write 
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Then,  the  volume  fraction  equation  for  fluid  (1)  can  be  rewritten  as 
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Equation  (182)  can  be  simplified  into  the  form 
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The  equation  for  can  be  derived  in  a  similar  manner. 


For  the  case  of  n  fluids,  the  k*  volume  fraction  equation  can  easily  be  deduced  from  the 


above  equation  and  is  given  by. 
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Improvement  #5;  Bounding  the  Volume  Fractions 

While  the  above  techniques  can  be  used  to  solve  the  volume  fraction  equations,  it  does  not 
guarantee  that  the  volume  fraction  values  are  bounded  (i.e.  between  0  and  1).  This  is  a  feature 
of  iterative  methods,  which  are  known  to  return  intermediate  values  that  violate  the  set 
bounds.  While  these  restrictions  can  be  explicitly  enforced  after  obtaining  the  solution  from 
the  discretized  equations,  a  number  of  techniques  were  developed  that  lead  to  the  implicit 
enforcing  of  these  constraints. 


The  procedure  developed  by  Carver  [74]  is  based  on  a  modification  of  the  under-relaxation 
procedure  of  Patankar  [16].  In  the  standard  under-relaxation  procedure  a  certain  proportion 
of  the  original  value  is  retained  in  the  solution.  Instead  of  solving  the  volume  fraction 


r 


A  Unified  Formulation  of  the  Segregated  Class  of Algorithms  for  Multi-fluid  Fluid  Flow  at  All  Speeds 


54 


equation  directly  for  r®  to  yield  under-relaxation  is  used  to  yield  a  value  j.(k)'mp 

+(1-P)  r^'‘^°,  where  P  is  between  [0,1]  and  is  the  solution  from  the  previous  iteration. 
Relaxation  may  be  imposed  after  the  solution  but  this  may  detract  from  the  integrity  of  the 
solution,  and  it  is  the  norm  to  pre-relax  the  equation,  making  the  relaxation  method  implied  in 
the  solution.  Equation  (25)  is  thus  rewritten  as: 


Yp  — 

P  NB 


^  A,n> 


(184) 


In  the  Carver  procedure,  the  value  of  over  each  control  volume  is  monitored  by 
explicitly  calculating  an  intermediate  value  r^^^'"*as: 


^  (A)int  NB 
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P 

If  rp*^'"‘>(l-8)  then  the  under-relaxation  factor  is  modified  to  the  form  p=MAX(l-rj*^,Y). 
The  parameters  s  and  y  are  small;  values  of  0.05  and  10'’ °  are  suggested  by  Carver.  The 
system  of  equations  for  each  volume  fraction  is  then  solved  implicitly  using,  for  every  control 
volume,  the  individually  assigned  relaxation  parameter. 


Solving  the  energy  equations 

The  solutions  of  the  energy  equations  follow  that  of  the  general  multi-fluid  scalar  equation. 
As  such,  nothing  new  needs  to  be  added  in  that  regard  (though  in  many  cases  coupling  of  the 
energy  equation  with  the  momentum  and  continuity  equations  is  beneficial). 
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Part  II  -  Geometric  Conservation  Based  Algorithm 


An  example  of  a  GCBA  is  the  original  IPSA,  developed  by  Spalding  [75],  which  introduces  a 
stronger  coupling  between  the  pressure  and  the  volume  fractions.  The  sequence  of  events  in 
the  Geometric  Conservation  Based  Algorithm  (GCBA)  is  as  follows: 

•  Solve  the  individual  mass  conservation  equations  for  volume  fractions. 

•  Solve  the  momentum  equations  for  velocities. 

•  Solve  the  pressure  correction  equation. 

•  Correct  velocity,  volume  fraction,  density,  and  pressure  fields. 

•  Solve  the  individual  energy  equations. 

•  Return  to  the  first  step  and  repeat  until  convergence. 

As  in  the  MCBA,  the  GCBA  uses  the  momentum  equations  for  a  first  estimate  of  velocities. 
However,  the  volume  fractions  are  calculated  without  enforcing  the  geometric  conservation 
equation.  Hence,  the  mass  conservation  equations  of  all  fluids  are  used  to  calculate  the 
volume  fractions.  The  pressure  correction  equation  is  based  on  the  geometric  conservation 
equation  and  is  used  to  restore  the  imbalance  of  volume  fractions.  The  errors  in  the 
calculated  volume  fractions  are  expressed  in  terms  of  pressure  correction  (P'  ,  which  is  also 
used  to  adjust  the  velocity  and  density  fields. 

Solving  for  volume  fractions 

The  solution  starts  by  estimating  the  volume  fractions  using  the  continuity  equations  given 
as: 


for  k  =  l..n  phases 


(187) 
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There  is  nothing  new  to  be  added  different  than  what  was  said  earlier  concerning  the  solution 
procedures  used  to  solve  the  above  set  of  equations  and  the  previously  described  techniques 
are  applicable. 

Solving  for  velocities 

Another  important  step  in  the  solution  cycle  is  to  obtain  an  initial  estimate  for  the  velocity 
fields  that  satisfy  the  momentum  equations.  For  that  purpose,  the  following  set  of  momentum 
equations  (given  here  for  completeness  of  presentation)  are  solved: 

uf  =  Up  ]-  VP  + 

m-  all  fluids^  k 

Again,  the  same  procedures  used  with  the  MCBA  can  be  applied  here  including  the  PEA  and 
SINCE  techniques.  Since  these  were  explained  earlier,  they  are  not  repeated  here. 

Solving  for  Pressure  Correction 

After  solving  the  continuity  equations  for  the  volume  fraction  fields  and  the  momentum 
equations  for  the  velocity  fields,  the  next  step  is  to  correct  the  various  fields  such  that  the 
volume  fraction  fields  satisfy  the  compatibility  equations  and  the  velocity  and  pressure  fields 
satisfy  the  continuity  equations.  For  that  purpose,  an  approach  similar  to  the  one  used  with 
the  MCBA  is  adopted.  The  difference  between  the  two  approaches  lies  in  the  constraint 
equation  employed  in  deriving  the  pressure  or  pressure  correction  equation.  In  the  MCBA, 
the  overall  mass  conservation  equation  was  utilized.  In  the  GCBA,  the  pressure  correction 
equation  is  derived  from  the  geometric  conservation  equation. 

To  start  the  derivation,  it  is  noticed  that  initially  the  volume  fraction  field  denoted  by  , 
does  not  satisfy  the  compatibility  equation  and  a  discrepancy  exists  i.e. 
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RESGp  =  \-Y,r^^^*  (188) 

k 

A  change  to  is  sought  that  would  restore  the  balance.  The  corrected  r  value,  denoted  by 

^(*)  _  ^(*)*  ^  fW  ^  jg  gm;}^ 

k  k 

or 


)=  1  - 

k  k 

Correction  to  the  volume  fraction,  ,  will  be  associated  with  a  correction  to  the  velocity, 
density,  and  pressure  fields,  ,  and  F  respectively.  Thus,  the  corrected  fields  are 

given  as: 
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Recalling  that  the  corrected  continuity  equation  of  phase  (k)  can  be  written  in  discretized 


form  as 
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a 


Its  expanded  form  becomes: 

(rr^rV'\pf°^pf^-{rypfT 

St 


(192) 


(193) 


or 


r 
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+  r^pV  +  rf'' pf")- 


a 


■  p  + 


rr/*)  u<*>* + u'*'*  1 

I  I  ^(4)*  . 


(194) 


Api:  + 


.S  !  = 


J 


Further  expansion  and  rearrangement  of  the  above  equation  yields 


a 


-#;vf  Qp=  - 


a 


(195) 


.r;"' p>''' u' ' 

Neglecting  second  and  third  order  terms  (i.e.  ,  and  rf"*  p^p^  ) 

the  above  equation  reduces  to: 


V  f’-  y'  f’'  +  A,[(r'‘'V<‘>u'*>.S  +  +  /’‘[/'‘V'*’)] 


vOW 


(196) 


-  jiy;v«>'n, = -- 


^  Jr,  p,  )  {r,  p,  )  _  Aj(/‘>V‘’V‘’')l+-A35.‘Vy’'n, 

Writing  as  a  function  of  F ,  similar  to  what  is  usually  done  in  a  simple-like  algorithm, 


the  correction  momentum  equations  become 

=  HP[u^*^']  -  -  r^*^'D^*VF 

Substituting  and  rearranging,  one  gets 


(197) 


_  (rr^,)^„.  __  ]  -  r<‘>'D“VF  -r'*''D'*Vr)s  +  r'^C/^V*'  (198) 


\0/d 


_  [fp  Pp  )  \rp  Pp  ) . .  ^itcy^jrjjjr^+ 

St 

Upon  discretizing  the  left-hand  side  of  the  above  equation,  it  can  be  rewritten  as 
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Ayf-;’  + 


NB  =  E,W,N,S,T,B 


■  Pp  Ap  r  p 


(199) 


{rrpf)-{rrdF 


a,- A,[(r'‘V‘’V‘’’)]+ 


By  noticing  that  the  left-hand  side  of  Eq.  (199)  resembles  the  volume  fi-action  equation  (Eq. 
39),  it  can  be  written  symbolically  as 


rf 


O‘>‘o  r  fHPru‘*’']-r‘*’‘D“’v/>''i  Ti 


d,  l_^(*)'D(*)VP' 


kVpfyirfp'^1 


where  ^  =  1  /  ylp*^ . 


Neglecting  the  correction  to  neighboring  cells,  equation  (200)  reduces  to: 


(200) 


r<‘V‘’’  .S  +  r<‘>W‘’'  I 

.<‘'=-4-1  ^  "  U-D-VP'  J  Jj 

^+V>  Pp  )  Yp  P’>  )  J 

Substitution  of  this  equation  into  the  geometric  conservation  equation  yields 


(201) 


f  r 


r^pf+A, 


arj^^VpA  Til 
.Sul 


SI -4  I 

T  (.■‘>vr)-(4py>r„  .  ■ 


!>  =  i?£5Gp 


Qp+A,[(r<V''°[/*>*)] 


(202) 


Writing  density  correction  in  terms  of  pressure  correction  (i.e.  ’  -  C^’F ),  equation  (202) 


is  transformed  to 
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I  'P 


I  I  ^  F,  +  1 1 

11^''  1 1 

i  I  •-  -'I  I 

{yW*  frW  I 

I +.VL_.^.j_VL.^_j_.  Q,+  A,[f  I 


(203) 


As  detailed  next,  the  above  equation  can  be  expanded  using  any  simple-like  algorithm  to 
yield  a  new  family  of  multi-fluid  flow  algorithms  (GCBA-SIMPLE,  GCBA-SIMPLEC, 
GCBA-PISO,...). 


The  GCB A  following  SIMPLE  (GCBA-SIMPLE):  Symbolic  Form 


Predictor: 


r;*>*  =  /f 


(204) 


Corrector: 


(205) 


(206) 


=  HP,[u<*>"]  -  =  HP,[u(^>*  +  -h  rf  ))WV,(P°  +  P  (207) 


(208) 


fp  —  JXp 


+  A,  [(r<^  )S  +  p^'^y  ]] 


Condition: 


j;f®'}=«£SG, 


(209) 
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I  I 

I  ,  dt  •-  -* 


II 


X  -Rp\  -  r^'^'D^*VP')s]j }  =  RESGp 

y[p  pp  )  yfp  pp  j . 


a 


frr^ 

St 


i^.+ApHr  p  ^ 

1,  Ap  [r  (hP[u^'’'  ]  -  Vp)s] 


Approximation: 

Neglect:  HP[u‘*’'], 

Approximate  Equation: 


(210) 


(211) 


(212) 


E« 


rp 


(k) 

^ '  * 


V 


f  ^ 

*  I 

A  Global  GCBA-SIMPLE  Iteration 


/  w»  w»\_  /  a)^wf ^1 
^  St  ^  A^[(r<^V*>°t/<^>*)]j|  -i?PAGp 


(213) 


•  Solve  implicitly  for  the  volume  fraction  fields. 

•  Solve  implicitly  for  u®  using  the  old  pressure,  density,  and  volume  fraction  fields. 

•  Calculate  the  D®  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCBA  following  SIMPLEC  (GCBA-SIMPLEC):  Symbolic  Form 


Predictor: 


(214) 
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(215) 


Corrector: 


(«“>', p"'.r 


(ky 


(216) 


+  F  (217) 

uf  =  HPp[u^*^'  ]-  r;*^*D^'VpP' - (218) 


Subtracting  HP/,[l]u^^^  from  both  sides,  one  gets 


u^'  -  HP^[1K^  = 


HP^[u(^>']-  HP^pjuf’  -P/^DfVpP' 


(219) 


1-  HP^[l]lu<;^'  =  mp[n^'‘^'-uf']-r]!‘^^DfWpP'-rj,‘‘^'D^^^Vpr-r^/^'D^PVpF  (220) 


rHPp[u<*^'  -u)f^'] 


Up  —  I 


1-HP,[1]) 


-n 


ar 


D 


(k) 


1-H'p,[1]) 


WpP' 


1 


(ky 

'p 


D' 


P) 


L 


(i-HP,[ii) 


VpF  -Vp 


iky 


jy 


ik) 


1-HP,[1] 


VpF 


J 


(221) 


f 


(*).  ^  HP^[u^^>'  nf]  _  jkrfyik)^  pf  _  ^<*)  d^V^F  -  D^^V^F 


\Up  = 


l-HP, 


.[1]) 


p  MJp  y  pr  fp  Ajp  y  pi  ip  M^p  y  pi 


i 


LE _ ^ 


^(ky  -  _  d(^)  ££- 

rp  -  Kp\  ^ 


pW'  +  [(r<^  u<*>')S  +  F*)*[/^*  I 


(222) 


Condition: 


{■«>'}= /!£SG, 


(223) 
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St 


r  r 


(k)*  (ky 

r  p 


(l-HP,[l]) 

[r,  p,  )  p,  )  ^  [(r<‘V*’'l^‘’*)] 


'll 

|l 

1 1|! 

.s!ll-=A£5Gp  (224) 


a 


1^  a 


11 


.  1^  V  llfJ— n,+  Aj(r 

2(4*'  a  f  av  ■  'j  |j._jfraG, 

[a,  P'V*’’  (hP[u<‘>'  -  u;'  ]  -  r'‘>'D<‘> VF')s] 


Approximation: 


(225) 


Neglect:  HP[u^'^' 

u(;>'  =  -r;*>*D<;V,P' 

Approximate  Equation: 


V 


St 


f  ^ 
*  i 


{r^pryir^^f) 

St 


Old 


.  Clp  +  A^[(r(*V*^V<'>‘)|  \-RESG, 


(226) 


(227) 


A  Global  GCBA-SIMPLEC  Iteration 


•  Solve  implicitly  for  the  volume  fraction  fields. 

•  Solve  implicitly  for  u®,  using  the  old  pressure,  density,  and  volume  fraction  fields. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  u®,  P,  and  p^''\ 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


r 
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The  GCBA  following  PRIME  (GCBA-PRIME):  Symbolic  Form 


Predictor: 


Corrector: 


(228) 

(229) 


(*)' 


u(i)**  =  ii(*)*  +  u(*)'  ^  -P°  P, 

=pik>  =  rW  + 


(230) 


.  =  HP,[u«**]  -  + p  (231) 

f 

u^'  =  HP^[u<*^*  -  +  nVpin^”^' ]  -  pP'  -  -  rfUfVpP' 

(232) 


Condition: 

Xfr}=^£SG, 

k 


(233) 


+ 


\  \  UL  ■- 

•  y  -«<*•!  A  r  <..  ,.rrHP["“>--u«>-]  +  HP[u<‘>-]^  ^1 

"k\  ^1  ^  J  j 

L  J' 


(234) 


^  /J  }.-PP5Gp 

^Ap[r‘ +  HP[u^*^']  -  r‘'^'D^*VP')s]Jj 


(235) 


A  Unified  Formulation  of  the  Segregated  Class  of  Algorithms  for  Multi-fluid  Fluid  Flow  at  All  Speeds 


65 


Approximation: 

Neglect: 

=>  pP' 

Approximate  Equation: 


a 


\  ^ 
*  l 


(rrpr)-{4^^pf) 

a 


Did 


■  +  A 


A  Global  GCBA-PRIME  Iteration 


(236) 


(237) 


•  Solve  explicitly  for  the  volume  fraction  fields. 

•  Solve  explicitly  for  u®,  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and  p®. 

•  Solve  Implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCBA  following  SIMPLEST  (GCBA-SIMPLEST):  Symbolic  Form 

Predictor: 


u*'  =  HPj'  [u“’']+  HP/[n“’']-  r“’'D™V,i>' 
Corrector: 


„(*)“  =  ,  P*  =  P’  +  F,  ^ 

p(k)*  -p(ky  j^pf^kY  ykr*  _  ^  ^.k)- ^ 


=  HP/’  [u<*>**]+  HP/[u<*^**  ]-  r«“DWV;,P* 

=  HP/  [u^*^*]+  HP/’[u^*^  ]+  HP/[u^*^’  ]+  HP/[u^*^  ] 
-  r/*^*D/V;,P°  -  /p^^D^p^pF  -  r/*''D/V^P°  - 


(238) 

(239) 


(240) 


(241) 
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f 


=HPp[u^*^']+ 

,w_^(k)]y  ^242) 


fliy  _  _  pw 

p  —  iVp 


Condition: 


Vrn. 

^  St 


pW' +  Ap[(r<  + 


•I 


(243) 


'  '  ^  ^  ^^i^  +  A^[rW*i7W’C</>F]+ 


St 


A, 


r  ...  ,,jHP[uW']  +  HPquW‘-u(^>]^ 


_^w*DWvP'  - 


.Si  I  !►  =  RESG, 


(244) 


a 


frrnA‘’ 


V  St 

^  ^(rrpn-(r^''pff 


p;  +  App^'t/^^^’Cj^^P']- App>V*>°(r^*^*D''VP')s||  = 


a 


Q,+  A,[(/V’'c/*0]+  I 


I 


(245) 


l^Ap  P*V*^°(hP^[**^*^*  -  +  HP[V*^']  -  r^*^'D^*VP')s]|  J 

Approximation: 


\-RESGp 


Neglect:  HP''[u^*>*- HP[u<*^'],  r^*>’D^*VP' 

=>  nf  =  -rf^^TifVpP' 

Approximate  Equation: 


(246) 


t  V 


-  ^<j  Rf  1^^^  ^  ^  +  A^  [(r ^  -  PP5G^ 


(247) 
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A  Global  GCBA-SIMPLEST  Iteration 

•  Solve  fort®  treating  implicitly  and  if  explicitly. 

•  Solve  for  u®  treating  implicitly  and  HP^  explicitly. 

•  Calculate  the  D^‘‘^  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and  p®. 

•  Solve  Implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCBA  following  SIMPLER  (GCBA-SIMPLER):  Symbolic  Form 


First  Predictor: 


MY  _ 


Calculate  the  coefficients  of  the  momentum  equations. 

First  Corrector: 

(uW'  p'  ) 

^  "1“  —  p^Y  _|_  j 

:.  =  HP^[u^'^‘  ]  -  VpP*  =  HPp[u^*^  +  (p°  +  P' 

=  HPp  ]- 

{ 


(248) 


(249) 

(250) 

(251) 


(252) 


kW'  =-pw 


'ifISh 

K  St 


Condition: 

Y,{f^'}=RESG, 


(253) 


f  (  P^*o 

I  I  ^^P^p 


2  \  -4*^1  Ap  -  r^'^^'D^*VP')s]j  \  =  RESGp 


(254) 


\Old 


+ 


<.>■,  P,  J  r,  Pf  )  A,[(r'‘V‘’’C^‘’-)] 


a 


A  Unified  Formulation  of the  Segregated  Class  of Algorithms  for  Multi-fluid  Fluid  Flow  at  All  Speeds 


68 


5t 

'  ( 


>  a,  +  A,[(/-'‘V'‘’V*'-)>  I 


a 


A,p‘>><*>‘  (hP[u''>'  ]  -  r'‘>'D<*' V/')S] 

Approximation: 

Neglect:  HP[u^'^^'], 

=> 

Approximate  Equation: 


\-RESGp 


V 


St 


FpFEp ]-  Ap [r Vf )s]  [  = 


f  " 
*  I 


{rrpV°)-{rrpf) 

St 


Did 


■  +  Ap 


'l-RESGf 

J 


Apply  correction  to  pressure,  density,  and  volume  fraction  fields. 
Second  Predictor: 


u 


(*)• 


Second  Corrector: 


(kT 


F'  =  F  F  F', 

=pW*  +pW,r<*)***  =r<*)** 


=  HP^  y  y^f^p  (P*  +  F' 

\ 

nf'  =  HP,[u'*^']  -  r^^DfVpF'  -  -  r^'^’DfVpP" 


D  “  “  r 


(255) 


(256) 


(257) 


(258) 


(259) 


(260) 


(261) 


Condition: 


Xfr'}=«£^G,  (262) 
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a 


^1 

II 

il! 


X ^ -4']  A, - r'*''D^^VP")s]j  ^  =  RESG, 

y-,  p,  )  v>  Pf  )  +A^|-(r<fl*y*)-(/*'')] 


(263) 


a 


f  r 


.Old 


/„W* 


(264) 


}-RESGp 


Approximation: 


Neglect:  HP[u^*^''], 

^nf'=-rr*T>f^pP"  (265) 

Approximate  Equation: 


a 


f  ^■ 

I  V 

Apply  correction  to  velocity  fields. 


f  f  /  (*)**  (ic)*\  /  (*)  aI 


(266) 


A  Global  GCBA-SIMPLER  Iteration 
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•  Solve  implicitly  for 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation  and  update  the  pressure,  density,  and  volume  fraction  fields. 

•  Solve  implicitly  for  using  the  new  pressure,  density,  and  volume  fraction  fields 

•  Solve  the  pressure  correction  equation  using  the  new  velocity  fieids  to  obtain  a  new 
pressure  correction  fieid. 

•  Correct  u® 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fieids. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCBA  following  PISO  (GCBA-PISO):  Symbolic  Form 


First  Predictor: 


First  Corrector: 


(267) 

(268) 


V)**  =uW*+uW',P*  =P“  +  F,  ' 

=p(k>  +  pW\p^y*  =  rW*  + 


(269) 


+  F  (270) 

f 


V  St 


)S  +  p^'^y  | 


Condition: 


I  I 


11^  L  J 

X  \  -4*^1  A^  -  r^*^'D^*VP')s]j  \  =  RESGp 


^1 

II 

.1! 


\Old 


[fp  Pp  )  yrp  Pp  )  ^  v*>*)] 


a 


(271) 


(272) 


(273) 
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(A)|  'P 

St 


-  y ^  <>!  ^  ^  ^  ^ 


+ 


^  Ap  [r  (hP[u<*^'  ]  -  r  Vf)s] 


\-RESGp 


Approximation: 

Neglect:  HP[u^*^'], 

u‘p^'  =  -r^*^*D^p  VpP' 

Approximate  Equation: 


St 

\  ^ 


m 


*  I 


/Mr  Mr\  /ik)  m'fi'k 
\p  Pp  }  \p  Pp  ) 


St 


•  Op  + 


.11 


\-RESGp 

J 


Second  Predictor: 


uT*  =  HP;[u^*^"]-rr*  Dlf^*‘VpP 
Second  Corrector: 


w 


u 


(ky 


****  =  =  P*  +  P", 


"  \ 


^p(k)-  -  p(kr  +  p(ky  ^^(k)****  _  ^(i)**»  +  jj.k)  j 


(274) 


(275) 


(276) 


.•.  = Hp;  [u^*>****  ]  -  Vp  (p*  +  p- 

u^***  =  HP;[u'*>“]  -  r;*>*’*D^*^“V^P‘ 

fu^'  =  HP;‘[uW”**  -  uW*‘ ]-  (r;*)’**  +  r;*)'))W**V^(P*  +P")+  rW“DW“V^P* 
j=  (hP;*[u<*>***  -u(*)**  +  u(*>"]-rW'iy;)"V^P"-rW**DW**V/,P"  -r(*)"DW"V^P* 
•’•  1  =  67*)p" 


(277) 

(278) 

(279) 

(280) 
(281) 


(282) 


I 


rW  _  vik)** 
r  p  ^Vp 


ll 


V  St 


pf  +  A^  [(p^>“  )s  +  r  J 


Condition: 

Efr'}=i-Z''r**=^‘5Gp 


(283) 
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p;'+  Jf 


Y  _/?(*)-!  A  ^  (*)“•  (*)-f HP**[u<*^“*  -  u<*>“  +  U<*>"]  )  1 


(284) 


11\R\ 


P;+  Ap  [r<*>***  1 


A^p)“y*)‘(/*)*«D(*)*VP")s] 


/ HP"[u<*)***  -  uW“  +  u(*>"]')  1 


.^(i)"D(*)**VP'' 


(285) 


-RESG, 


Approximation: 


Neglect;  HP”[u^*^“* -u^*^** +  u^*^’’], 
=  -ri*"D«"V,P" 


Approximate  Equation: 


f  P*)***Q  r 

,A^p^>**V*>*(r(^>***D(*>‘*VP")s]  J  J 


f  f  (j.(ky**(ky  \_  (  ik)  (i)Y>« 

\P  Pp  )  Vp  Pp  ) 


) -  JiFSG^ 


(286) 


(287) 


A  Global  GCBA-PISO  Iteration 


A  Unified  Formulation  of  the  Segregated  Class  of Algorithms  for  Multi-fluid  Fluid  Flow  at  All  Speeds 


73 


•  Solve  implicitly  for 

•  Solve  implicitly  for  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and 

•  Solve  implicitly  the  energy  equation  and  update  the  density  fields. 

•  Solve  explicitly  for 

•  Solve  the  momentum  equations  explicitly  and  calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and 

•  Return  to  step  one  and  iterate  until  convergence _ 


The  GCBA  following  SIMPLEX  (GCBA-SIMPLEX):  Symbolic  Form 
Predictor: 


(288) 

(289) 


Corrector: 


('  W  r>i  (*) 

u  ,P  ,p  ,r 


,  ,  P*  =  r  +  F, 


^pik>  ^pW^^ikY-  ^  j.(kT  +  j.(kYj 


(290) 

-  rW”D?>V^P*  =  HP^[u<*)*  +  -  (rf *  +  +  F  (291) 


f 

((*)*/-» 


Condition: 

k 


(292) 


(293) 


p  p 


St 


>]|  =  RESGj, 

fo‘vy)-fo‘v?’r 


(294) 


St 
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Approximation: 

Neglect:  D^p\pF  ,  and  let 

u^p^'  =  HPp[u^*^'  ]-  pP' 

Assume  that  the  pressure  difference  local  to  the  velocity  is  representative  of  all  pressure 
differences  i.e.  HPp[-/'^*D<*^^-^VP']  =  ,  thus: 


(295) 

(296) 


-r^*^*Djf^^^(V;,P')=  -(VpP')HPp[/*’*D‘*^^'^]  -  (VpF) 


=  HPp[P^^  +  r^lip 

Approximate  Equation: 


Jkr-r.(k)sx-,  ,  „(i)*rk(i) 


(297) 

(298) 


p  p 


5t 


f  ^ 
*  I 


St 


Did 


Qp  +  A^ 


(299) 


\-RESGp 
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A  Global  GCBA-SIMPLEX  Iteration 


•  Solve  implicitly  for 

•  Solve  implicitly  for  u®,  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  implicitly  for  the  fields. 

(k)SX 

•  Solve  the  pressure  correction  equation  using  these  D  fields. 

•  Correct  u®  P,  and 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  GCBA  following  SIMPLEM  (GCBA-SIMPLEM):  Symbolic  Form 

First  Predictor: 


Calculate  the  coefficients  of  the  momentum  equations. 

First  Corrector: 


,  P*  =  P°  +P',  '' 


(300) 


(301) 
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+^r’)Dp^V/,(P°  +  P' 


f 


u^p^'  =  HPp[u‘^’']-  rj'>’D^'VpP'-r^*^'D^*VpP°  -  VpP' 


.-.  ^  =  C[^^F 


V  St 


I'll 


pf  +  Ap  [(r<^ 


'll 


f  C<*> 

^  Uwl  ^^P^P 


St 


,.(*X  -_vW 
'p  ~  P^p 

Condition: 

k 

I  I  />  +  App*>V^*^Cf  P']+ 

I  I  II 

.-.  X  ^  -P^*'|  Ap  P>>^*^XhP[«'*"]  -  VP'  -  VP')s]|  ^  =  RESG, 

'  Qp+Ap[pV^^°t/^>°)]  J 

I  V'"  Pp"  )  f  Cl  +A  I 

Ap  (hP[u^*^’  ]  -  D‘^  VP')s]  J  J 

Approximation: 

Neglect:  HP[u^*>'],  r<*)'D<*VP' 

^  uf  =  D^p  VpP' 

Approximate  Equation: 

X{<(^^^^^^i?  +  A,p>'l7<*V/’/'']-A,[/‘'y*’‘(r<‘’'D'‘>V/>')s||  = 

f  f  rP^y^- A'‘'>'f'‘‘  ^1 

_  Vp  Pp  )  ^p  Pp  )  n,  +  A  ^  -  RESG, 

*  [  V  ^ 


'  (i 
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(304) 
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Second  Predictor: 

pP  (310) 

Second  Corrector: 


No  corrector  stage. 

A  Global  MCBA-SIMPLEM  Iteration 


•  Solve  implicitly  for  ~ 

•  Calculate  the  fields  based  on  values  from  the  previous  iteration. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,  and 

•  Calculate  new  and  fields. 

•  Solve  implicitly  for  0^“^^  using  the  new  fields. 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

» Return  to  the  first  step  and  iterate  until  convergence. _ 

As  for  the  MCBA,  most  of  the  GCBA  algorithms  have  not  been  extended  and  tested  in  multi¬ 
fluid  flow  situations.  As  such,  the  only  way  to  discuss  their  merit  is  by  assuming  that  they 
behave  as  in  single-fluid  flows.  Since  such  a  discussion  was  given  earlier,  it  will  not  be 
repeated  here. 


The  Expanded  Form  of  the  Pressure-Correction  Equation 


The  expanded  form  of  the  pressure  correction  equation,  applicable  to  all  algorithms,  will  be 
presented.  For  that  purpose,  let  r^’'^  and  denote  values  from  the  previous  iteration  or 
from  a  previous  corrector  step,  then,  the  pressure  correction  equation  becomes 


5t 


p;+  A^pV'^(/V>V(p'))s||  = 


nphase  I  (  w- 

_  Vp  Pp  )  Pp  ) 

*=1  [  V  7  J 

In  this  form,  the  resemblance  between  this  equation  and  the  MCBA  pressure  correction 
equation  is  obvious.  The  discretization  of  the  above  equation  yields 


(311) 


a;  Fp  =  ApP^ +a;;f^+ a;;p; + Afp^ + a^p;  +  4  +  p, 


(312) 


where 
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A^J  =  Ap  +  AC  +  AC  +  Af  +  Ar'  +  A^'  + 


f  ik)  J^(k)r\ 

f  rp  Kp  Qp  t/f  ^  ^ 


{ (^{k)  (k)  _  (k)old  (k)old^ 

I  \rp  Pp  Tp  Pp 


[  +  dpfu]!^^)  J 


\-RESGp 


(313) 


Following  the  calculation  of  the  pressure  correction  field,  Up^ ,  dp^ ,  and  are  obtained 


using  the  following  equations 


uW'  =  -r<*>DWVp(P') 

.•.^W'=C^*)F  (314) 

Improvement  #6 

The  pressure  correction  equation  was  derived  in  the  previous  section  by  neglecting  the 
HpY'‘^  term.  Following  the  SIMPLEC  methodology,  a  better  approximation  can  be 
achieved  by  adding  and  subtracting  Hp\i\rf  fi’om  the  left-hand  side  of  equation  (199), 
which  results  in  neglecting  a  smaller  term  {HpY''^  -rf  \  ■  With  this  approximation,  the 


pressure  correction  equation  becomes; 
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(1-H,[l])rr  = 


- 
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+r 


ay 


(r,  p,  ;-  (r,  p,  )  ^  Aj(r'‘>V‘>V*>‘)]+ jy;v;‘>'a,  ) 


or 


(315) 


Cmr,  rr<*>y«'(HP[u<*>']-r<‘>'D<*>VF  - r<*>'D<*> VF')s'll 
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+rWu^>‘r  w 


J| 


\Old 


Where 


(316) 


D(k) 

p(k)  _  ^ 


(317) 


The  pressure  correction  equation  is  obtained  by  substituting  Eq.  (316)  into  the  geometric 
conservation  equation.  The  expanded  form  of  the  resulting  pressure  correction  equation  is 

easily  obtained  from  Eqs.  (312)  and  (313)  by  simply  substituting  Rp^  for  RpK  Furthermore, 
depending  on  the  approximation  made  to  the  term,  a  new  family  of  GCBA,  similar 

to  the  one  detailed  above,  can  be  obtained. 


Part  III  -  Comparing  GCBA  and  MCBA  Formulations 


Scaling  GCBA  to  Single-Fluid  Flow  Simulations 

Since  MCBA  algorithms  are  derived  through  direct  extension  of  the  single  fluid  algorithms, 
they  can  scale  down  automatically  to  handle  one-fluid  simulations.  On  the  other  hand, 
because  GCBA  algorithms  are  based  on  geometric  conservation  (as  opposed  to  mass 
conservation)  their  scaling  to  one-fluid  simulations  is  not  obvious.  Such  a  property  is  useful 
in  the  sense  that  coding  for  single  and  multi-fluid  models  would  follow  the  same  structure. 
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To  show  how  multi-fluid  GCBA  algorithms  scale  down  to  single-fluid  flow  simulations, 
equation  (201)  is  reproduced  as  follows: 


I 


L 


.S  +  rW*C/(*>y*>' 


(318) 


\OW 


+ 


A(*)* 


J 


For  the  case  of  one  fluid  flow  r^^^=l=Cte,  thus  rp^  removing  from  equation  (318) 
and  setting  to  1  yields: 


0  =  -R 


ik) 


-r<*^'D^*VF)S  + 


(319) 


Removing  ftP  and  rearraging,  the  pressure  correction  equation  transforms  to 


The  expansion  of  the  above  equation  is  straightforward  and  leads  to  the  pressure-correction 
equation  of  a  single-fluid  flow. 

Relating  GCBA  and  MCBA 

While  the  derivations  of  GCBA  and  MCBA  are  based  on  different  paradigms,  it  is  shown  in 
this  section  that  the  GCBA  formulation  leads  to  a  weighted  pressure  correction  equation  that 
has  close  similarity  with  the  MCBA  formulation  to  which  improvement  #2  has  been  applied. 

The  respective  GCBA  and  MCBA  (with  improvement  #2)  pressure  correction  equations,  as 
derived  earlier,  are  as  follows: 
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(*)  'p  '^p  ^‘‘P 

''  I  5t 

nphase  f  ^ 
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ir,  p,  ;  p,  )  ^  A,[(/‘V<‘'C/'‘>)] 


1 

\-RESGp 
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(321) 


f  [  Sk>f<t)rj(k)p,^  \  (k)  1] 

vJ  “  +  A  ^  A  '  w^/;.(*rpWVP>^  ■ 


(322) 


Upong  comparing  the  two  equations  it  is  clear  that  /?p  *  plays  the  role  of  the  reference 
density,  i.e.  a  weighing  faetor  for  the  respective  phasic  continuity  equation.  As  such,  the 
GCBA  pressure  correction  equation  is  very  similar  to  a  weighted  MCBA  pressure  correction 
equation.  The  weighing  procedure  (normalization)  is  done  automatieally  based  on  the  local 
strength  of  the  inflow  to  the  control  volume  (sinee  =  \/Ap^of  the  volume  fraction 
equation)  thus,  as  one  gets  «  Rf^ .  This  treatment  yields  a  more  robust 

behaviour  since  large  fluid  density  differences  will  not  mean  that  conservation  for  the  lighter 
fluid  is  lost  due  to  numerical  errors. 


The  above  equations  also  differ  slightly  in  the  source  term,  where  an  additional  entry 
(RESGp)  is  included  in  equation  (321)  to  aecount  for  residuals  of  the  geometrie  conservation 
equation.  Hence,  the  family  of  MCBA  algorithms  can  be  viewed  as  a  subset  of  the  GCBA 
family,  and  could  be  recovered  by  setting  the  volume  fraction  corrections  to  zero  (r^*^  =  0). 
As  such,  codes  based  on  the  MCBA  can  easily  cater  for  GCBA  and  vice  versa. 

The  advantage  of  the  GCBA  over  the  MCBA  is  in  its  attempt  to  correct  both  the  velocity  and 
volume  fraction  fields.  However  there  are  still  limitations  pertaining  to  the  ealculation  of  the 
volume  fraction  field,  namely  that  the  ealculation  is  done  without  aceounting  for  the 
influence  of  the  volumes  fraetions  of  the  different  fluids  on  each  other.  Thus  while  the 
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pressure  change  will  tend  to  balance  the  volume  fractions,  it  will  not  do  so  in  an  optimal  way. 
In  addition,  all  volume  fractions  will  change  in  the  same  direction  in  response  to  pressure 
adjustment,  which  is  not  always  logical. 

Closing  Remarks 

The  segregated  class  of  single-fluid  flow  algorithms  was  extended  to  predict  multi-fluid  flow 
at  all  speeds.  The  formulation  was  done  using  a  unified,  compact,  and  easy  to  understand 
notation.  Depending  on  the  constraint  equation  used  to  derive  the  pressure  or  pressure 
correction  equation,  the  extended  algorithms  were  shown  to  fall  under  two  categories  that 
were  denoted  by  the  Mass  Conservation  Based  Algorithms  (MCBA)  and  the  Geometric 
Conservation  Based  Algorithms  (GCBA).  The  differences  and  similarities  between  the  two 
categories  were  explained.  In  addition,  several  techniques  developed  to  promote  and 
accelerate  the  convergence  of  these  algorithms  were  also  discussed. 
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